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Abstract. We present the results of MHD simulations of low mass protoplanets interacting with turbulent, 
magnetised protostellar disks. We calculate the orbital evolution of 'planetesimals' and protoplanets with masses 
in the range < m p < 30 M®. The disk models are cylindrical models with toroidal net-flux magnetic fields, 
having aspect ratio H/r = 0.07 and effective viscous stress parameter a ~ 5 x 10 -3 . 

A significant result is that the m p — 'planetesimals', and protoplanets of all masses considered, undergo stochas- 
tic migration due to gravitational interaction with turbulent density fluctuations in the disk. For simulation run 
times currently feasible (covering between 100 - 150 planet orbits), the stochastic migration dominates over type I 
migration for many models. Fourier analysis of the torques experienced by protoplanets indicates that the torque 
fluctuations contain components with significant power whose time scales of variation are similar to the simulation 
run times. These long term torque fluctuations in part explain the dominance of stochastic torques in the models, 
and may provide a powerful means of counteracting the effects of type I migration acting on some planets in 
turbulent disks. The effect of superposing type I migration torques appropriate for laminar disks on the stochastic 
torques was examined. This analysis predicts that a greater degree of inward migration should occur than was 
observed in the MHD simulations. This may be a first hint that type 1 torques are modified in a turbulent disk, 
but the results are not conclusive on this matter. 

The turbulence is found to be a significant source of eccentricity driving, with the 'planetesimals' attaining 
eccentricities in the range 0.02 < e < 0.14 during the simulations. The eccentricity evolution of the protoplanets 
shows strong dependence on the protoplanet mass. Protoplanets with mass m v = 1 Mq attained eccentricities in 
the range 0.02 < e < 0.08. Those with m v = 10 Mj reached 0.02 < e < 0.03. This trend is in basic agreement 
with a model in which eccentricity growth arises because of turbulent forcing, and eccentricity damping occurs 
through interaction with disk material at coorbital Lindblad resonances. 

These results are significant for the theory of planet formation. Stochastic migration may provide a means of 
preventing at least some planetary cores from migrating into the central star due to type I migration before they 
become gas giants. The growth of planetary cores may be enhanced by preventing isolation during planetesimal 
accretion. The excitation of eccentricity by the turbulence, however, may act to reduce the growth rates of 
planetary cores during the runaway and oligarchic growth stages, and may cause collisions between planetesimals 
to be destructive rather than accumulative. 
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1. Introduction 

The continuing discovery of extrasolar planets by the ra- 
dial velocity and transit techniques has generated renewed 
interest in the theory of planet formation (e.g. Mayor & 
Queloz 1995; Marcy, Cochran, & Mayor 2000; Vogt et al. 
2002; Santos et al. 2003). In the so-called core instability 
model of planet formation, gas giant planets form through 
the build-up of a rocky and icy core of ~ 15 Earth masses, 
which then undergoes gas accretion resulting in a giant 



planet (e.g. Bodenheimer & Pollack 1986; Pollack et al. 
1996). The gas accretion stage of this process is normally 
believed to take a few million years, during which time 
the planetary cores slowly accrete gas through quasi-static 
settling, until the planet mass is around 50 Earth masses. 
Beyond this mass the gas accretion rate increases dramat- 
ically and gas is able to accrete onto the planet through a 
circumplanetary disk at the rate supplied by the protostel- 
lar disk (Papaloizou & Nelson 2005). An alternative model 
for giant planet formation is that gravitational fragmenta- 
tion of a fairly massive protostellar disk can form Jovian 
mass planets directly (e.g. Boss 2001). This model, how- 
ever, is unable to account for terrestrial planets, and is 
unlikely to explain the existence of Uranus and Neptune. 
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Since giant planets form in a gas-rich protoplanetary 
disk, the gravitational interaction between planets and 
disks can play an important evolutionary role. This in- 
teraction has been the subject of many studies over the 
last couple of decades. In the usual picture, a protoplanet 
exerts torques on a protostellar disk through the exci- 
tation of spiral density waves at Lindblad resonances, 
and possibly through interaction at corotation resonance 
(e.g. Goldreich & Tremaine 1979, 1980; Lin & Papaloizou 
1979; Papaloizou & Lin 1984; Ward 1986, 1997; Tanaka, 
Takeuchi & Ward 2002). The spiral waves carry an as- 
sociated angular momentum flux, which is deposited in 
the disk material where the waves damp, leading to ex- 
change of angular momentum between protoplanet and 
disk. The disk lying exterior to the protoplanet orbit ex- 
erts a negative torque on the planet, and the interior disk 
exerts a positive torque. The negative outer disk torque 
usually dominates and the protoplanet migrates inward. 
This is referred to as type I migration. For protoplan- 
ets of ~ 15 Earth masses, the migration time is between 
10 4 and 10 5 yr (Tanaka, Takeuchi & Ward 2002), which 
is much shorter than the estimated gas accretion phase 
of giant planet formation (Pollack et al. 1996). Taken at 
face value, this presents a serious problem for the core- 
instability model of gas giant planet formation. In a recent 
paper, Papaloizou & Nelson (2005) examined the possi- 
bility of shortening the formation time scale of gas giant 
planets by reducing the dust opacity in the outer radia- 
tive zone of the forming planets. Although the formation 
time can be shortened significantly, they concluded that 
the type I migration time is always shorter than the planet 
formation time. This analysis, however, pertains only to 
smooth, laminar disk models. 

Until quite recently most models of viscous accretion 
disks used the Shakura & Sunyaev (1973) a model for the 
anomalous disk viscosity. This assumes that the viscous 
stress is proportional to the thermal pressure in the disk, 
without specifying the origin of the viscous stress. Balbus 
& Hawley (1991) showed that significant angular momen- 
tum transport in weakly magnetised disks could arise from 
the magnetorotational instability (MRI). Non linear nu- 
merical simulations performed using a local shearing box 
formalism (e.g. Hawley & Balbus 1991; Hawley, Gammie, 
& Balbus 1996; Brandenburg et al. 1996) confirmed this 
and showed that the saturated non linear outcome of the 
MRI is MHD turbulence having an effective viscous stress 
parameter a of between ~5x 10 -3 and ~ 0.1, depending 
on the magnetic field configuration. Global simulations of 
MHD turbulent disks (e.g. Armitage 1998; Hawley 2000; 
Hawley 2001; Steinacker & Papaloizou 2002; Papaloizou 
& Nelson 2003) confirm the picture provided by the local 
shearing box simulations. 

An obvious question is how the interaction of form- 
ing planets with protoplanetary disks changes when the 
disks are turbulent. Papaloizou & Nelson (2003 - here- 
after PN2003) examined and characterised the turbulence 
obtained in a variety of MHD disk models. Nelson & 
Papaloizou (2003 - hereafter NP2003) examined the in- 



teraction between a global cylindrical disk model and 
a massive (5 Jupiter masses) protoplanet. A similar 
study was undertaken by Winters, Balbus, & Hawley 
(2003). Papaloizou, Nelson, & Sncllgrove (2004 - here- 
after PNS2004) performed global cylindrical disk simu- 
lations and local shearing box simulations of turbulent 
disks interacting with protoplanets of different mass. The 
main focus of that paper was to characterise the changes 
in flow morphology and disk structure as a function of 
planet mass, and to examine the transition from linear to 
non linear interaction leading to gap formation. Nelson & 
Papaloizou (2004 - hereafter NP2004) examined the mi- 
gration of low and high mass protoplanets in turbulent 
disks by examining the time evolution of the torques ex- 
erted on the planets by the disks. They noticed that low 
mass planets experienced strongly varying torques due to 
interaction with the turbulent density fluctuations, and 
suggested that such planets would undergo stochastic mi- 
gration rather than monotonic inward migration normally 
associated with type I migration in laminar disks. The 
planets, however, were maintained on fixed circular orbits 
and so their migration could not be examined directly. 

In this paper we present the results of simulations of 
low mass protoplanets embedded in turbulent, magnetised 
disks, and allow the planet orbits to evolve due to inter- 
action with the disk. In order to sample the range of out- 
comes, we consider a small ensemble of planets in each 
simulation. We find that in all simulations performed, the 
forces experienced by the protoplanet are highly variable, 
and as a result the protoplanet orbits evolve similarly to a 
random walk. A simple analysis suggests that the heavier 
planets we consider ought to undergo inward migration 
predominantly due to the underlying type I torques, but 
the simulations on the whole are not in agreement with 
this prediction. Fourier analysis of the torques experienced 
by the planets indicates that they experience stochastic 
forces with a broad range of associated time scales of vari- 
ation, ranging from the planet orbital period to the sim- 
ulation run time itself. These long time scale variations 
clearly contribute to the long term stochastic behaviour 
of the planets observed in the simulations. There is some 
evidence that the underlying type I migration is modified 
in turbulent disks, but this is not conclusive. 

The results also indicate that MHD turbulence is able 
to drive significant growth of eccentricities for low mass 
objects. In particular, planctesimals and planets with 
masses similar to the Earth can obtain eccentricities e ~ 
0.1. This clearly has potentially important consequences 
for planetary accumulation models. Heavier planets attain 
lower eccentricities apparently because their interaction 
with the disk at coorbital Lindblad resonances causes ec- 
centricity damping (e.g. Artymowicz 1993; Papaloizou & 
Larwood 2000). 

The plan of the paper is as follows. In section 2 we de- 
scribe the governing equations. In sections 3 and 4 we de- 
scribe the numerical method and the system of units used. 
The initial and boundary conditions used in the simula- 
tions are described in section 5, and the simulation results 
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are presented in section 6. We discuss the simulation re- 
sults in sections 7, 8, and 9, focusing in particular on the 
issues of the balance between type I and stochastic migra- 
tion, and the eccentricity evolution of planets and plan- 
etesimals. We summarise the paper and draw conclusions 
in section 10. 



2. Equations of Motion 

The equations of ideal MHD written in a frame rotating 
with uniform angular velocity Q/z, with z being the unit 
vector along the rotation axis assumed to be in the vertical 
direction, are the continuity equation 



dp 
dt 



+ V • pv = 0, 



the equation of motion 



9v „ Vp „^ (V x B) x B 
— + v • Vv + 2fi f zxv = - V$ + -r 2 (2 

at p Airp 



and the induction equation 
(9B 

-=Vx(vxB) 



(3) 



where v, P, p, and B denote the fluid velocity, pressure, 
density and magnetic field respectively. The potential $ = 
®rot + ®g contains contributions due to gravity, and 
the centrifugal potential $ rot = -(l/2)0^|z x r| 2 . 

The gravitational potential has contributions from a 
central mass M* and ./V planets with masses m P i. Thus in 
cylindrical coordinates (r,<f),z), with the planets located 
at (r pi , 4> p i,0) and the star located at the origin of the 
coordinate system, the gravitational potential is <J>g = 



4\: 



and 



Gm r , 



\Jr 2 + r 2 pl - 2rr pl cos(</> - (j) pi ) + b 2 



(4) 



(5) 



Here, as in the papers PN2003, NP2003, PNS2004 and 
NP2004, we have neglected the dependence of the gravita- 
tional potentials on z along with the vertical stratification 
of the disk, for reasons of computational speed. Thus the 
simulations are of cylindrical disks (e.g. Armitage 1998, 
2001; Hawley 2000, 2001; Steinacker & Papaloizou 2002). 
To model the effects of the reduction of the planet poten- 
tial with vertical height, we have incorporated a softening 
length b in the potential. 

We use a locally isothermal equation of state in the form 



P = c s {rf 



(6) 



with c s denoting the sound speed which is specified as a 
fixed function of the cylindrical radius r. Expressions (1) 
- (6) give the basic equations used for the simulations. 



2.1. Planet Orbital Evolution 

In this work we allow the planet orbits to evolve due to 
the gravitational forces they experience from the disk as 
well as the central star. For a single simulation we consider 
the evolution of either six or three planets concurrently, 
located at different positions within the disk. These plan- 
ets do not interact with each other gravitationally, but 
affect each others evolution indirectly through perturba- 
tions they make to the disk structure. In the case of low 
mass planets, the effects of the turbulence will be of con- 
siderably greater significance. 
The equation of motion for planet i is: 



(1) dt 



GNL 



r P i +Ui - 20/z x v p , +0/z x (O/z x r pi ). (7) 



The acceleration due to the disk is given by 



Jv 



p(r){r pi -r) 



(r 2 + r 2 pi - 2rr pi cos{<j> - (j) pi ) + b 2 ) 3 / 



;dV (8) 



where the integral is performed over the disk volume. Note 
that the planet potentials are also cylindrical in that there 
is no vertical component of gravity. During the simulations 
we monitor the torque per unit mass experienced by the 
planets, which is defined by 



^pi X fdi • 



(9) 



Note that the torque per unit mass is independent of the 
planet mass, and so may be calculated for zero-mass 'plan- 
etesimals' and for finite mass planets. 

The calculations are performed in a uniformly rotat- 
ing frame with angular velocity f2 / . The planet equations 
of motion are evolved using a simple leap-frog integra- 
tor, with correct centering employed for inclusion of the 
Corioli's force. 

3. Numerical Method 

The numerical scheme that we employ is based on a spa- 
tially second-order accurate method that computes the 
advection using the monotonic transport algorithm (Van 
Leer 1977). The MHD section of the code uses the method 
of characteristics constrained transport (MOCCT) as out- 
lined in Hawley & Stone (1995) and implemented in the 
ZEUS code. The code has been developed from a version 
of NIRVANA written originally by U. Ziegler (Ziegler & 
Yorke 1997). 

4. Computational Units 

We use units in which the central mass = 1, the 
gravitational constant G = 1, and the radius R = 1 
corresponds to the radial location of the inner bound- 
ary of the computational domain. The unit of time is 
fl^ 1 = y/GM^/R 3 , although we report our results in units 
of the orbital period at the disk inner edge, which we de- 
note as P(l) = 2iril~ 1 . Planets are positioned in the disk 
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Fig. 1. This figure shows the time evolution of the volume 
averaged total a value in the disk. Note that the time is ex- 
pressed in units of the orbital period at the inner edge of the 
disk. 

between radii r V i — 2.2 - 3.2. We note that the orbital 
period at r pl = 2.2 is P(2.2) = 3.26P(1) and at r pi = 3.2 
we have P(3.2) = 5.72P(1). Simulation runs times are 
typically ~ 500 P(l), corresponding to ~ 153P(2.2) and 
~ 90P(3.2). If we take the computational radius r = 2.5 to 
correspond to 5 AU, then a simulation run time of 500P(1) 
corresponds to ~ 1335 yr, for a solar mass central star. 

5. Initial and Boundary Conditions 

The simulations presented here all use the same under- 
lying turbulent disk model. This model has a constant 
aspect ratio H/r = c(r)/(rO) = 0.07, where f2 is the disk 
angular velocity measured in the inertial frame. The in- 
ner radial boundary of the computational domain is at 
Ri n = 1 and the outer boundary is at R ou t — 5. The 
simulations were performed in the rotating frame with 
flf = 0.28668, which is the Keplerian angular velocity 
at a radius of r = 2.3. 

The boundary conditions employed are very similar 
to those described in PN2003. Regions of the disk in the 
vicinity of the inner and outer boundaries were given non- 
Keplerian angular velocity profiles (uniform in the simu- 
lations described here) that are stable to the MRI, and 
which have large values of the density in order to main- 
tain radial hydrostatic equilibrium. These regions act as 
buffer zones that prevent the penetration of magnetic field 
to the radial boundaries, thus maintaining the initial value 
of net magnetic flux in the computational domain. The in- 
ner buffer zone runs from r = 1.2 to r = R m . The outer 
buffer zone runs from r = 4.5 to r = R ou t- 

The creation of a relaxed turbulent disk model in which 
planets could be immersed to examine their orbital evolu- 
tion was constructed using a multi-stage procedure. This 
is because the initial relaxation of a turbulent disk model 
can significantly modify the initial surface density struc- 
ture, due to radial and temporal variations in the magnetic 



Fig. 2. This figure shows the evolution of the radial variation 
in the total a value in the disk. The solid line corresponds to 
a time of 15.6 orbits at the inner disk edge, the dashed line 
corresponds to 31.4 orbits, the dotted line to 71.2 orbits, and 
the dash-dotted line to 188.1 orbits. 

stresses during the relaxation. It is desirable to minimise 
this effect so that a reasonably smooth profile is used for 
the disk in which the planets are inserted. Also, we wish 
to obtain a turbulent disk model in which the volume av- 
eraged a value is ~ 5 x 10 -3 , as expected for protostellar 
disks, and the magnetic field strength in the disk was ad- 
justed to achieve this. 

The calculations presented in NP2004 used an initial 
setup where Ri n — 1 and R ou t — 8. The numerical resolu- 
tion used was greater than that used here ([N r , Nj,, N z ]= 
[450, 1092, 44], as opposed to [264, 608, 44] used in this 
paper), but those simulations could only be run for ~ 25 
planet orbits. The higher resolution used in NP2004 meant 
that zero net flux magnetic fields could be used, requiring 
a self-sustaining dynamo to maintain the field. The lower 
resolution adopted for the simulations in this paper leads 
to the requirement that net flux magnetic fields be used 
to sustain the turbulence over long time periods. This is 
because the lower resolution simulations cannot maintain 
an active turbulent state for zero net flux fields. The use of 
net flux fields guarantees the continued existence of mag- 
netic field in the simulation domain, which helps maintain 
the turbulence. A modest magnetic flux is used to obtain 
a disk with modest turbulent stresses. 

In order to obtain the desired disk model, the following 
procedure was adopted: 

(i). The disk was initiated with a toroidal magnetic field in 
a finite annulus in the disk between the radii 1.5 < r < 4. 
The initial field varied as B^r) = Bo/r, and Bo was de- 
fined such that the initial ratio of volume integrated mag- 
netic pressure to thermal pressure was l/(/3) = 2.73 x 10 -3 
(where j3 = 8ttP/B 2 ). The initial density distribution in 
the disk was p(r) = po/r away from the buffer zones, with 
po defined such that the disk model would contain ap- 
proximately 0.06 M Q interior to 40 AU (assuming r = 2.5 
corresponds to 5 AU). In other words, the basic model is 
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chosen to be around three times more massive than the 
minimum mass solar nebula model. The disk model was 
relaxed for a period of 191.7 orbits at the disk inner bound- 
ary. 

(it). The disk model was restarted at this point with the 
initial density distribution re-established, but with all 
other quantities remaining the same after the disk was 
relaxed for 191.7 orbits. The model was then run for a 
further 58.3 orbits. 

(Hi). The model was restarted with the initial density 
distribution re-established, and with the magnetic field 
strength being reduced throughout the disk by a factor of 
1/V2. The model was run for a further 37.2 orbits. 
(iv). The model was restarted after a further reduction 
in the magnetic field strength by a factor of l/y/2, and 
with the initial density profile having been re-established. 
The model was relaxed for a further 32 orbits. The fi- 
nal model has a volume averaged (a) = 5.04 x 1CT 3 , and 
l/(/3) = 0.012, which are similar to those obtained in the 
models presented in PN2003, PNS2004 and NP2004. 

The disk model has a full 2ir azimuthal domain, and 
vertical domain from z m i n = —0.14 to z rnax = 0.14. 
Periodic boundary conditions were employed in the verti- 
cal and azimuthal directions. 

Gravitational softening of the protoplanet was em- 
ployed in all simulations, with the softening parameter 
b = 0.33-ff and H being evaluated at each planet position. 

6. Simulation Results 

6.1. Turbulent stresses 

In order to describe average properties of the turbulent 
models, we use quantities that are both vertically and az- 
imuthally averaged over the (<f), z) domain (e.g. Hawley 
2000). The vertical and azimuthal average of Q is defined 
through 



5v<j,{r, <fi, z, t) = v<t,(r, <j), z, t) - v<j,(r, t). 



(15) 



Q(r,t) = 



J pQdzd(f> 
J pdzdip 



(10) 



The average is taken over the full 2ir in azimuth. The disk 
surface density is given by 



2tt/ 



pdzdcf). 



(11) 



The vertically and azimuthally averaged Maxwell and 
Reynolds stresses are defined as follows: 



T M (r, t) = JhrE ( B r {r^z,t)B,{r^z,t) ^ 



and 



T Re (r, t) = 27r£(5u r (r, 0, z, t)5v<f,(r, 4>, z, t). 



(12) 



(13) 



Here the velocity fluctuations Sv r and Sv^ are defined 
through, 



The Shakura & Sunyaev (1973) a stress parameter appro- 
priate to the total stress is defined by 



a(r, t) 



Trs — Tm 



2ttE(P/p)' 



(16) 



Sv r (r, <j), z, t) = v r (r, <f), z, t) — v r (r, t), 



(14) 



6.2. Disk Model 

The method used to set up the turbulent disk model is 
described in section 5. The subsequent evolution of the 
volume averaged a value is shown in figure 1, for the du- 
ration of the simulation performed with zero-mass planets 
described below. Values of the volume averaged stress pa- 
rameter (a) ~ 5 x 10~ 3 are sustained throughout the du- 
ration of the simulation, which is approximately 500 orbits 
at the disk inner edge. Snapshots of the radial distribution 
of a are shown in figure 2, showing that a magnetically 
active disk is sustained, at least for radii beyond r > 2. 

Figure 3 shows snapshots of the midplane density for 
the model with six 10 M ffi protoplanets embedded. These 
images show that the density perturbations induced by 
planets of this mass and lower are exceeded by those gen- 
erated by the turbulence. The consequences of this on 
the planet orbital evolution are described below for planet 
masses ranging between < m p < 30 M®. The effects of 
the turbulence on protoplanet orbital evolution are essen- 
tially stochastic, and for this reason the simulations were 
performed with six or three planets in order to sample the 
range of possible outcomes. 

6.3. Zero mass protoplanets - 'planetesimals' 

The orbital evolution of the six 'zero-mass planets' is 
shown in figure 4. These planets began with orbital radii 
r p i — 2.2, 2.4, 2.6, 2.8, 3.0, 3.2 and randomised azimuthal 
positions. From now on we shall refer to these zero-mass 
objects as planetesimals, as they represent the evolution 
of bodies whose gravity is too weak to substantially per- 
turb the disk structure, and for which no underlying type I 
migration occurs. We note that planetesimals experience 
gas drag, which is able to modify their orbits. We have 
neglected this effect here, but models of planetesimal evo- 
lution including gas drag will be presented in a future 
paper. 

The left hand panel of figure 4 shows the evolution 
of the planetesimals' semimajor axes, with each line rep- 
resenting the evolution of a different planetesimal. It is 
clear from this plot that the gravitational interaction of 
the planetesimals with the turbulent density fluctuations 
causes the semimajor axes to change stochastically. This 
results in the planetesimals undergoing a 'random walk' in 
orbital location, with variations of the semimajor axis at 
the 5-10 percent level being observed over time periods 
corresponding to <~ 100 planetesimal orbits. 

The right hand panel shows the evolution of the plan- 
etesimal eccentricities. Not surprisingly the effect of the 
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Fig. 3. This figure shows snapshot images of the disk midplane density for the run with m p — 10 planets. It is clear that 
the turbulent density fluctuations are typically larger than the spiral wakes generated by the planets in this case. The times 
corresponding to each image are shown at the top right of each panel, in units of the orbital period at the disk inner edge. 



turbulence is to excite significant orbital eccentricity in 
the planetesimal population. Typical values obtained are 
e ~ 0.03, but with values ranging from 0.01 < e < 0.13 
over the run times of ~ 100 planetesimal orbits. 

The results presented in figures 4 suggest that MHD 
turbulence will have a significant effect on the orbital evo- 
lution of planctcsimals in protoplanetary disks, which will 
in turn have important repercussions for models of planet 
formation. In particular, the effects of stochastic migra- 
tion will help to prevent the isolation of forming planetary 
embryos, and reduce the effects of orbital repulsion and 
planetesimal shepherding during oligarchic growth (Ida & 
Makino 1993; Kokubo & Ida 1998; Thommes et al 2003). 
This may allow the growth of planetary cores to proceed 
to the 10 - 15 M e range, required for the formation of 
giant planets via the core instability scenario (e.g. Pollack 
et al. 1996). However, it is also clear that MHD turbulence 
may hinder the growth of planetary cores by exciting sig- 
nificant orbital eccentricities of the planctcsimals, reduc- 



ing the effects of gravitational focusing. For planetesimals 
in the 1-10 km range, it is unlikely that gas drag will 
be effective in reducing these eccentricities because of the 
dominant contribution from the turbulence as a source of 
excitation. Under these conditions an alternative source 
of eccentricity damping may be required, such as inelastic 
collisions between planetesimals, leading also to fragmen- 
tation and more effective gas drag for the smaller frag- 
ments. Of course, it then remains an open question how 
planetesimals can form in the first place in such a turbu- 
lent environment. This may be related to the existence of 
'dead zones' in which the ionisation fraction is too low to 
support MHD-driven turbulence (e.g. Gammie 1996). 

6.4. Finte mass planets 

The simulations presented in the next five sections (6.5 - 
6.9) were performed for either six or three planets placed 
in turbulent disks. These planets do not interact with each 
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Semirnajor axes versus time Eccentricity versus, time 
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Fig. 4. This figure shows the evolution of the planet semimajor axes (left panel) and eccentricities (right panel) for the run 
with planet mass set to zero. It is clear from this plot that the time varying gravitational field of the disk causes stochastic 
migration of the 'planetesimals', and significant eccentricity growth. 



other gravitationally, as the purpose of these simulations 
is to sample the effects of turbulence on single planets. The 
planets may, however, affect each other indirectly through 
their influence on the disk. For planets whose orbits are 
well separated, this influence appears to be insignificant 
compared to the turbulence. In some simulations, how- 
ever, there are planets that approach one another and 
occassionally cross orbits, increasing the mutual indirect 
influence. Such simulations will no longer be sampling the 
effects of MHD turbulence on isolated planets. We note 
that this issue is likely to be most important for the more 
massive planets. 

The planets in the simulation involving the 30 M ffi ob- 
jects do not approach each other. The simulation including 
the 10 Mq protoplanets does involve close encounters. We 
have examined the torques experienced by the planets, 
along with their semimajor axis and eccentricity evolu- 
tion, and can detect no obvious sign that the evolution is 
being modified by the close encounters. We are confident 
that the main results of this paper are not contaminated 
by this effect. 

6.5. m p = 1 Earth mass protoplanets 

Figure 5 shows the orbital evolution of the 1 Mq proto- 
planets. The dotted lines in the left panel show the evo- 
lution of planet orbits for simulations in which six 1 M e 
planets were evolved in a laminar disk model that was 
otherwise identical to the turbulent disk model. The solid 
lines, corresponding to planets in the turbulent disk, show 
similar stochastic migration as observed for the planetes- 



imals. Indeed the protoplanet initially located at r p — 2.4 
migrates inward for ~ 70 orbits measured at the inner 
disk edge at first, before migrating outward for a sustained 
period of ~ 130 orbits, modifying its semimajor axis by 
almost 20 percent during this time. The other planets are 
also seen to migrate in a quasi-random fashion, but over 
less extreme distances. 

The evolution of the eccentricities is plotted in the 
right hand panel. A tendency for eccentricity growth to 
significant values for some of the protoplanets is observed, 
in line with the results for the planetesimals. It is inter- 
esting to observe however, that the peak eccentricities ob- 
tained are somewhat less than those obtained in the plan- 
etesimal cases, suggesting that the gravitational interac- 
tion with the disk is providing eccentricity damping to 
counteract the effects of the turbulent excitation. This is 
a trend that continues as the protoplanet mass increases, 
as described in the following sections. 

6.6. m p = 3 Earth mass protoplanets 

The orbital evolution of the 3 protoplanets is illus- 
trated in figure 6. The left panel shows a similar evolu- 
tion for the semimajor axes as has been described for the 
planetesimals and 1 M ffi protoplanets. The semimajor axes 
are again observed to undergo modification in a stochastic 
manner, with variations at the 5 percent level being ob- 
served. The dotted lines show the results of a simulation 
in which six 3 M ffi protoplanets were evolved in a laminar 
disk. The monotonic inward migration of each of these 
planets is clearly discernible, suggesting that the planets 
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Fig. 5. This figure shows the evolution of the planet semirnajor axes (left panel) and eccentricities (right panel) for the run 
with m p = 1 M®. The dotted lines show the evolution of planet orbits for simulations in which six 1 M® planets were evolved 
in a laminar disk model that was otherwise identical to the turbulent disk model. The migration rates obtained are in quite 
good agreement with the rates predicted by Papaloizou & Larwood (2000) for planets with softened potentials migrating in 
laminar disks. The solid lines represent the planets in the turbulent disk, and show that the planets undergo migration similar 
to a random walk for the duration of the simulation, with no clear tendency for the planets to migrate inward or outward. 



do not greatly influence each others interaction with the 
disk. Inward migration in this case occurs close to the 
expected type I rate. The simulated torques are smaller 
than those presented in Papaloizou & Larwood (2000) by 
a factor of 0.71. 

The eccentricity evolution is illustrated in the right 
panel of figure 6. Significant variations in the eccentrici- 
ties are again observed, with peak values of e ~ 0.05 being 
excited. It is interesting to compare the evolution with the 
planetesimal and 1 M e cases, as there is a clear trend to- 
ward obtaining lower eccentricities with increasing planet 
mass. It is likely that the eccentricity evolution is deter- 
mined by a balance between stochastic forcing due to the 
turbulence, and eccentricity damping due to interaction 
with the disk at coorbital Lindblad resonances (Goldreich 
& Tremaine 1980; Artymowicz 1993). This point is exam- 
ined and discussed further in section 7. 

6. 7. m p = 5 Earth mass protoplanets 

The orbital evolution of the six 5 M ffi protoplanets is 
shown in figure 7. The semirnajor axes in the left panel 
show a general inward trend, which might at first be 
thought to be evidence that the usual type I migration 
is becoming dominant over the stochastic forcing as the 
planet mass increases. Closer inspection of the figures 
shows that the protoplanets are still subject to substantial 
random forcing. Furthermore the 10 M ffl case discussed be- 
low suggests that the general inward motion of the planets 
in figure 7 is probably a statistical effect rather than ev- 



idence for type I migration overwhelming stochastic forc- 
ing. It is worth noting, however, that for all cases consid- 
ered in which the planets have masses m P i > 1 M®, the 
inner most planet migrates inward at a rate similar to or 
greater that the expected type I rate. This arises because 
the planet moves into a region of the disk where it is less 
turbulent due to the initial set up, and where the density 
increases during the course of the simulation due to disk 
accretion. 

The eccentricity evolution shown in the right panel of 
figure 7 continues to show the trend of lower eccentricity 
for higher mass planets. Here the peak value for the ec- 
centricity obtained is e ~ 0.04, but with eccentricities in 
the range 0.01 < e < 0.02 being more typical. 

6.8. m p = 10 Earth mass protoplanets 

Figure 8 shows the orbital evolution of the six 10 M ffi pro- 
toplanets. Examination of the left panel shows that the 
planets' orbital evolution remains dominated by stochas- 
tic forcing for this mass, with the semirnajor axes of three 
of the planets showing net increases for the duration of the 
simulation. This suggests that the general inward motion 
for the 5 M® planets seen in figure 7 is a statistical ef- 
fect. The dotted lines in the left panel of figure 8 show the 
evolution of six 10 M ffi planets in a laminar disk. The ex- 
pected inward migration is observed for each planet close 
to the predicted rate (Papaloizou & Larwood 2000), show- 
ing that the planets have little effect on each others orbital 
evolution due to their perturbation of the disk. 



R.P.Nelson: Low mass planets in turbulent disks 



9 



Seminncijor axes versus time Eccentricity versus, time 




100 200 300 400 500 100 200 300 400 500 

Time {Orbits) Time {Orbits) 

Fig. 6. This figure shows the evolution of the planet semimajor axes (left panel) and eccentricities (right panel) for the run with 
m p = 3 M ffi . The figures show that the planets undergo migration similar to a random walk for the duration of the simulation, 
with no clear tendency for the planets to migrate inward or outward. 
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Fig. 7. This figure shows the evolution of the planet semimajor axes (left panel) and eccentricities (right panel) for the run 
with m p = 5 M®. The figures show a general inward trend for the migration, but this appears to be a statistical fluke as the 
m p — 10 cases do not. Note that the scale on the y axis of the right panel has changed compared with that in similar 
figures 4-6. 



The eccentricity evolution is shown in the right panel 
of figure 8. It is clear that the eccentricities remain quite 
small in this case, with peak values reaching e ~ 0.025. 
More generally, however, the eccentricities remain at ~ 
0.01. The overall trend of eccentricity with planet mass 
remains such that considerably less eccentricity growth 
is observed for higher mass planets. The indication is 



that disk-planet interaction at coorbital Lindblad reso- 
nances causes eccentricity damping, and for planet masses 
of Tripi ~ 10 M ffi this dominates over turbulent forcing. 
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Fig. 8. This figure shows the evolution of the planet semimajor axes (left panel) and eccentricities (right panel) for the run with 
m v = 10 M®. The figures show that the planets undergo migration similar to a random walk for the duration of the simulation, 
with no clear tendency for the planets to migrate inward or outward. Note that the scale on the y axis of the right panel has 
changed compared with that in similar figures 4-7. 



6.9. m p = 30 Earth mass protoplanets 

We consider the orbital evolution of three 30 M ffi pro- 
toplanets, in contrast to the six planets considered for 
smaller mass objects, since they exert larger perturbations 
on the disk. Snapshots of the midplane density distribu- 
tion are shown in figure 9, which show that the planets 
are now massive enough to excite density waves of simi- 
lar amplitude to the background turbulent fluctuations. A 
similar trend was noted by NP2004. 

The semimajor axis evolution is shown in the left panel 
of figure 10. The inner most planet is observed to migrate 
inward, such that after 400 orbits it has reached the same 
location as predicted for a laminar disk. This is again re- 
lated to the fact that the inner disk is less turbulent. The 
outer two planets, however, show a significant deviation 
from the trajectories followed by planets in a laminar disk 
model. The middle planet undergoes noisy inward migra- 
tion, but at a rate substantially less than that obtained 
for the laminar model. The outermost planet undergoes 
essentially no net migration for the duration of the model 
run time. It is clear that stochastic effects due to the tur- 
bulence continue to have an impact on the migration of 
30 M e protoplanets. 

The eccentricity evolution is shown in the right panel of 
figure 10. The values obtained are similar to those for the 
10 M ffi planets, indicating little growth of the eccentricity 
because the damping induced by the underlying type I 
resonant disk interaction. 



7. Stochastic torques and type I migration 

We now examine the torques experienced by the proto- 
planets during the simulations in more detail. In view of 
the large number of planets considered, we restrict our 
discussion to a few specific examples which are illustra- 
tive of the range of behaviour observed in the simulations 
as a whole. We also discuss the implications for type I 
migration of low mass planets in turbulent disks. 

7.1. Time evolution of stochastic torques 

Figure 11 shows four examples of the time evolution of 
the torque per unit mass experienced by the planets in 
the simulations. The torque per unit mass is defined by 
equation (9), and is presented in the units described in sec- 
tion 4. Moving from left to right and from top to bottom, 
the panels show results from simulations with: (m P i = 0, 
r pi = 2.6); {m pl = 1, r pi = 2.4); (m pl = 3, r pi = 2.8); 
(m p i = 10, r p i = 2.6), where the radii, r p i, refer to the ini- 
tial orbital radii of the planets. In each panel, the torque 
on the planet due to the inner disk is shown by the black 
(blue) line, that due to the outer disk is shown by the light 
grey (green) line, and the dark grey (red) line represents 
the total torque. In all cases the torque is a highly vari- 
able quantity, as discussed previously in NP2004. On the 
same scale, the type I torque due to an equivalent lami- 
nar disk is ~ 1.5 x 10~ 6 (^j^j ■ This value was obtained 
from the simulations of low mass planets embedded in 
laminar disks represented by the dotted lines in figures 5 
- 8. Comparison with the migration time given by equa- 
tion (32) of Papaloizou & Larwood (2000) shows that our 



R.P.Nelson: Low mass planets in turbulent disks 



11 



30 Earth mass planets 



30 Earth mass planets 




6- >, □ 



-2 




-4 



r7 



30 Earth mass planets 



30 Earth mass planets 




-4 



-2 



j i_ 

4 




-4 



-2 



Fig. 9. This figure shows snapshot images of the disk midplane density for the run with m p — 30 planets. It is clear that 
the turbulent density fluctuations are of similar amplitude to the spiral wakes generated by the planets in this case. The times 
corresponding to each image are shown at the top right of each panel, in units of the orbital period at the disk inner edge. 



simulated torques are smaller by a factor of ~ 0.71 com- 
pared to their fit to linear calculations. Stochastic varia- 
tions relative to type I torques therefore typically range 
between 10 and 100 for planet masses in the range 1 to 
10 M®, but with peak fluctuations being up to four times 
larger than these values. As was noted in NP2004, these 
large fluctuations make it difficult to separate the torques 
due to the inner and outer disk when inspecting figure 11 
for planet masses m p < 10 M®. 

The time evolution of the running means of the torques 
displayed in figure 11 are shown in figure 12, along with 
the running means of the torques from the equivalent lam- 
inar disk simulations. In all cases, the running means do 
not converge to the equivalent type I laminar disk torques 
during the simulations, as one would expect from inspect- 
ing the migration histories presented in figures 4 - 8. If 
good convergence of the running means to the appropri- 
ate type I torque is possible, then these figures show that 
runs times considerably longer than those presented here 



will be required to achieve it. Such runs are currently not 
feasible computationally. 

This latter point concerning convergence is illuminated 
to some degree by figure 13, which shows the power spec- 
trum of the temporal evolution of the total torques plot- 
ted in figure 11. In these figures, a frequency of between 
0.2 - 0.3 represents the orbital periods of protoplanets lo- 
cated at orbital radii between r P i =2.23 - 2.92. The power 
spectra are obtained by Fourier analysis of the time evolu- 
tion of the total torques generated during the simulations. 
These torques are output every ten time steps during the 
simulations, with the time interval between each output 
being At. If N is the total number of data points then we 
can define discrete frequencies by 



2nn 
N At 



(17) 
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Fig. 10. This figure shows the evolution of the planet semimajor axes (left panel) and eccentricities (right panel) for the run 
with m p — 30 M ffi . The figures show that the outer two planets undergo migration that is substantially different to the laminar 
disk runs. The inner planet migrates primarily inward as the disk turbulence is reduced in the inner regions. Note that the scale 
on the y axis in the right panel differs from that of similar plots pertaining to lower mass protoplanets. 



where n — ~N/2, N/2. The discrete Fourier transform 
of the torques is then 

1 tf ~ 1 

H{u n ) = — ^ T k ex P (iu n tk) (18) 

where tk is the time corresponding to data point k, and 
Tk is the torque. We define an amplitude by 

A(u n ) = K)| 2 + l# (-^n)| 2 . (19) 

A(u>„) is plotted against |a> n | in figure 13. 

The power spectra show that there is significant signal 
in the long term variation of the torques, with the total 
run time corresponding to a frequency of 0.002 in these 
plots. This explains why the averaged torques in figure 12 
do not converge to a well defined value: the stochastic 
torques contain contributions with significant amplitude 
whose associated time scale of variation is similar to the 
simulation run time. 

The origin of these low frequency fluctuations is un- 
clear. One possibility is that global communication across 
the disk on the viscous evolution time may lead to mod- 
ification of the local turbulence on that time scale. This 
may in turn feed into the torque fluctuations experienced 
by embedded protoplanets. The disks we consider in this 
paper do not have a very large range in radii, but for the 
disk parameters H/r = 0.07 and a — 5 x 10 -3 the vis- 
cous communication time between the region where the 
planets are located (centred around r = 2.8) and the disk 
outer edge (r = 4.5) is r v ~ 8800 orbits. Given that the 
simulations have been run for 500 orbits, communication 



between neighbouring regions of the disk may be influ- 
encing the turbulence on this time scale. One test of this 
would be to compare the power spectra of fluctuations in 
global disk calculations with local simulations performed 
in the shearing box. If the local simulations do not show 
the long term fluctuations then this could be taken as 
evidence that global communication is important. If low 
frequency fluctuations are observed, then it would indi- 
cate that these fluctuations are the property of the locally 
operating dynamo. 

The torque frequency distributions corresponding to 
figure 11 are plotted in figure 14. These were calculated 
by sampling the torque experienced by the planets every 
10 time steps during the simulations. The torques were 
then binned with bin widths AT = 10~ 5 , and the num- 
ber in each bin counted. Although not strictly Gaussian, 
these distributions are similar to Gaussian profiles with 
standard deviation between <jt — 1 - 2 x 10~ 4 . 

7.2. Comparing stochastic migration and type I 
migration torques 

In a previous study, NP2004 considered the orbital evo- 
lution of low mass planets in turbulent disks, and treated 
the issue of torque convergence as a simple signal-to-noise 
problem. The basic assumption here is that the torque ex- 
perienced by the planet, T(t), consists of a linear superpo- 
sition of a rapidly varying term, Tf(t), caused by turbulent 
fluctuations, and a constantly acting type I torque < T >: 

T(t) =<T> +T f (t). (20) 
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Fig. 11. This figure shows the time evolution of the torques per unit mass for four different planets selected from four different 
simulations. Moving from left to right and from top to bottom the planet masses and initial orbital radii are: (m pi = 0, r pi = 2.6); 



(rripi = 1, r pi = 2.4); (m pi = 3, 



2.8); (m P i — 10, r V i — 2.6) The light grey (green) line corresponds to the torque exerted 



on the planet by the outer disk, the black (blue) line shows the torque due to the inner disk, and the dark grey (red) line shows 
the total torque. Moving from m v = to m p = 10 Me we observe that the inner and outer disk torques tend to separate, but 
that the torque fluctuations due to the turbulence remain significantly greater in amplitude than the mean torque. 



Assuming that the stochastic torques are Gaussian dis- 
tributed (the Central Limit Theorem tells us that even if 
the torques are not normally distributed, their cumulative 
effects will be similar to normally distributed torques), the 
time average of equation (20) becomes 

T=<T (21) 
yttot 

Here <jt is the standard deviation of the torque amplitude, 
and Uot is the total time elapsed, measured in units of 
the characteristic time for the torque amplitude to vary. 
Convergence of the torques toward the underlying type I 
value is expected to occur once the two terms on the right 
hand side become equal. 

Using a very simple by-eye inspection of 'torque ver- 
sus time' plots, similar to figure 11, NP2004 estimated 
the time for torque variation and the typical amplitude 
of the fluctuating torques. Using these in equation (21) 
they derived a time for torque convergence of ~ 70 planet 
orbits for a 10 M ffi protoplanet. Following a similar proce- 



dure, inspecting a blow-up of figure 12 leads to a by-eye 
estimate of the torque variation time ~ half a planet or- 
bital period (the same as used by NP2004). Inspecting the 
torque distrbutions in figure 14 gives a value of the torque 
standard deviation <jt — 1-5 x 10~ 4 . Inserting these val- 
ues into equation (21), and using a type I torque value 

of 1.5 x 10~ 6 (j^j, gives an estimated time for torque 

convergence equivalent to 50 planet orbital periods when 
the protoplanet mass is 10 M®. The modest difference in 
this estimate compared to NP2004 arises because the es- 
timates of the standard deviation of torque fluctuations 
differ. For planets distributed in radius as they are in our 
simulations, this translates into a convergence time of be- 
tween 163 - 286 orbits measured at r = 1. This crude 
analysis predicts that we should see inward migration for 
the 10 M ffi planets in figure 8, which is clearly not the 
case. 

Treating the issue of torque convergence using the 
above approach has a number of problems. The main one 
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Fig. 12. This figure shows the time evolution of the running mean torques per unit mass for the four simulations shown in 
figure 11. The planet masses and initial orbital radii are indicated at the top right of each panel. The light grey (green) line 
corresponds to the torque exerted on the planet by the outer disk, the black (blue) line shows the torque due to the inner 
disk, and the dark grey (red) line shows the total torque. Also plotted with black lines are the running means obtained from 
equivalent simulations with viscous, laminar disk models. Here the upper line is due to the inner disk, the lower line is due to 
ther outer disk, and the middle line is the total torque. 



is illustrated by figure 13, which clearly demonstrates that 
there is no single frequency associated with the torque 
variation, and furthermore indicates that a simple by-eye 
examination of torque versus time plots can be misleading 
about which time scale is dominant. A more sophisticated 
treatment should take account of the fact the stochas- 
tic torques have significant amplitude distributed across a 
broad range of frequencies. We have taken a different ap- 
proach to addressing the issue of torque convergence that 
ought to be more accurate. 

In order to test whether the evolution of low mass 
planets in turbulent disks can be described as a super- 
position of stochastic torques and constantly acting type 
I torques, we have performed a series of N-body simu- 
lations in which particle orbits were evolved under the 
influence of prescribed forces. These included stochastic 
forces, which were calculated using a Fourier analysis of 
the migration torques obtained in the MHD simulation 
with planet masses m p i = (i.e. the 'planetesimal' run 



shown in figure 4). This Fourier analysis allowed forces 
from the actual MHD simulations to be incorporated into 
simple integrations of particle orbits. The Fourier trans- 
form of the torques as a function of time is given by equa- 
tion (18) in section 7.1. The reconstructed torque at time 
t is then given by the inverse transform: 

N-l 

T(t) = exp {-i(J n t) (22) 

n=0 

This expression can be calculated for any value of t re- 
quired, and so provides a means of including the stochastic 
torques in an N-body simulation. We note, however, that 
the torques calculated during the MHD simulations were 
only output once every ten time steps, rather than at ev- 
ery time step. This means that our reconstructed torques 
are not exactly the same as the torques experienced by 
the planets during the MHD simulations. The variation of 
the torque experienced by a planet over ten time steps is 
not great (the torque shows significant variation over ~ 
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Fig. 13. This shows the power spectrum of the total torques plotted in figure 11. The main point to be taken from these plots 
is that significant power exists in frequencies that correspond to the longest time scales in the simulation, suggesting that low 
frequency variations in the torques experienced by the planets are important. 



one thousand time steps), so we expect the reconstructed 
torques to be a good approximation to the real torques. 

In addition to including stochastic force, type I torques 
taken from equation (32) of Papaloizou & Larwood (2000) 
were included, for which the inward migration rate is lin- 
early proportional to the planet mass. We note that com- 
paring torques from laminar disk simulations, and those 
given in Papaloizou & Larwood (2000), show that those 
obtained in the simulations are smaller by a factor of 
~ 0.71, so we multiplied the type I torques used in the 
N-body simulations by this factor. 

The particle orbits were integrated using a simple leap- 
frog integrator. The time step used was held constant, and 
was chosen to be equal to the mean time step size which 
arises from the Courant condition in the MHD simula- 
tions. The time step calculated during the MHD simula- 
tions varied by only a few percent throughout the simula- 
tions. 

The results of these runs are shown in figure 15, where 
the planet mass is indicated in the upper right corner of 
each panel. For zero-mass planets, orbit evolution similar 



to that shown in figure 4 was obtained. This indicates that 
the sampling of torques after every ten time steps in the 
MHD simulations leads to reconstructed torques that are 
quite accurate. For successively higher masses, these orbit 
trajectories were modified by having a small inward drift 
superposed. For planet masses m p i > 10 M®, a definite in- 
ward drift that overwhelms the stochastic migration can 
be observed. It is worth noting that the inward drift ob- 
served for the 10 M e cases in figure 15 is not matched by 
the actual MHD simulation presented in figure 8. 

We should point out that we have only sampled six 
realisations of the turbulence when performing the cal- 
culations presented in figure 15. On this basis alone it 
would be wrong to draw too strong a conclusion about 
the role of underlying type I torques in the MHD simula- 
tions. Inspection of figures 4-8, and figures 13 and 14 show 
that the torque distributions can differ substantially from 
run to run, such that the results shown in figure 15 may 
simply be a reflection of a minority of possible outcomes. 

To test this further we ought to perform additional 
simulations of MHD turbulence with zero-mass test par- 
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Fig. 14. This shows the frequency distribution of the torques plotted in figure 11. The plots show that the stochastic torques 
experienced by the planets have distributions that are fairly close to being Gaussian with standard deviation o(T) ~ 1— 2 
xlO" 4 . 



tides being used to sample the torque fluctuations. This, 
however, would be computationally expensive. The re- 
maining simulations we have presented, for finite mass 
planets, contain the effects of turbulent fluctuations and 
spiral density waves (not always visible against the turbu- 
lent backdrop) excited by the planets. As such they cannot 
be used to sample cleanly the role of turbulence-induced 
fluctuations, but including them in the current analysis is 
still useful. 

A 10 Mg planet in a laminar disk migrates a distance 
of A r ~ 0.15 in 500 orbits (as shown by the dotted lines in 
figure 8). We now consider the effect of superposing such 
a drift on the migration histories of all planets plotted in 
figures 4-8 whose orbits started between radii r V i = 2.6- 
3.2. Of this sample of twenty planets, a subset of eleven 
show net outward migration prior to the superposition. Of 
these eleven, only one would show net outward migration 
after superposition (this would be the 10 M® planet that 
was initiated at r p i — 3.0), two would show no net migra- 
tion, and eight would show a net inward drift. Based on 



this apparent tendency for the simulations as a whole to 
show inward migration when type I drift is superposed, 
which differs from the outcome observed in figure 8 for 
which the planet mass is 10 M®, there is a hint that the 
effects of type I migration may be diminished in a tur- 
bulent disk. This may arise because of density and pres- 
sure perturbations in the vicinity of the planet modifying 
the excitation of spiral density waves and the usual bias 
between inner and outer disk. At present, however, the 
results of the simulations are sufficiently ambiguous for 
such a conclusion to be premature, ft remains possible 
that the fluctuation spectrum associated with the orbital 
evolution shown in figure 8 is such that it is simply able to 
overcome the effects of type I migration for the duration 
of the simulations due to significant contributions from 
low frequency components. Resolution of these questions 
requires simulations to be run for longer duration, and for 
customised simulations to be performed that examine in 
detail the exchange of angular momentum between planet 
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Fig. 15. This figure shows the evolution of semirnajor axes from simulations in which the stochastic torques generated by the 
MHD simulation with m p = were reconstructed using Fourier analysis, and type I migration torques were included using a 
scaled version of equation (32) from Papaloizou & Larwood (2000). For planet masses lower than 10 M®, stochastic evolution 
of the orbits is dominant. For m p — 10 M®, the effects of the superposed type I migration torques are dominant. 



and turbulent disk. These calculations will be the subject 
of future publications. 

8. Eccentricity Evolution 

In order to examine the origin of the behaviour of the or- 
bital eccentricities observed in figures 4 - 8, we performed 
simulations similar to those described at the end of sec- 
tion 7.2. Particle orbits were evolved using a simple leap- 
frog integrator, and with prescribed forces included rep- 
resenting stochastic eccentricity forcing and eccentricity 
damping associated with interaction at coorbital Lindblad 
resonances. The stochastic forces were calculated using a 
Fourier reconstruction of the radial and azimuthal forces 
experienced by the 'planetesimals' in the MHD simula- 
tion for which the planet mass m p i = 0, results for which 
are shown in figure 4. The eccentricity damping was in- 
cluded by using equation (38) of Papaloizou & Larwood 



(2000). As already described, the source of this damping is 
primarily interaction with the disk at coorbital Lindblad 
resonances (e.g. Artymowicz 1993). 

We calculated orbital evolution for planets with masses 
m p i = 0, 1, 3, 5, 10, and 30 M e as in the MHD simula- 
tions. The masses simply control the degree of eccentric- 
ity damping, whereas the eccentricity driving is indepen- 
dent of mass. A selection of the resulting eccentricity be- 
haviours are shown in figure 16. These figures show that 
good agreement is obtained with the evolutionary trends 
found in the full MHD simulations, suggesting that the 
eccentricity evolution obtained in those runs was indeed 
due to stochastic forcing by the turbulence, and damp- 
ing due to resonant disk interaction. It is interesting to 
note that the agreement with the eccentricity evolution is 
better than that obtained when considering the superposi- 
tion of type I torques and stochastic forces in section 7.2. 
The reason for this is simply that the type I migration 
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Fig. 16. This figure shows the evolution of eccentricities from simulations in which the stochastic torques generated by the MHD 
simulation with m p — were reconstructed using Fourier analysis, and eccentricity damping was included using equation (38) 
from Papaloizou & Larwood (2000). The planet mass in each case is indicated in the top right of each panel. A clear trend in 
decreasing eccentricity may be observed as the masses of the planets increase, in broad agreement with the full MHD simulation 
results presented in sections 6.3-6.9. 



torques rely on there being an asymmetry between the 
inner and outer disk contributions, which may not oc- 
cur as the density field near the planet varies erratically. 
The eccentricity damping, however, is largely produced 
by coorbital Lindblad torques located near the orbit of 
the planet. These should always be present throughout 
the MHD simulations, and operate as a constant source 
of damping, although the rate of damping may vary with 
time due to fluctuations induced in the disk structure by 
the turbulence. 

9. Discussion 

9.1. Long term evolution of planetesimal orbits 

Here we consider only objects for which gravitational fluc- 
tuations due to the turbulent disk dominate the dynamics 
rather than gas-drag forces. Plantesimals with radius 1 



- 10 km fall into this category. The migration histories 
of the simulated planetesimals plotted in figure 4 show 
that such bodies will undergo substantial migration via a 
random walk if protoplanetary disks are globally turbu- 
lent. The mean migration distance of the planetesimals 
in figure 4 is Ar ~ 0.09 over a time of 500 orbits mea- 
sured at r = 1. If we make the simplifying assumption 
that the migration distance scales according to a random 
walk, so that Ar a \/t, where t is the time elapsed, and 
further take the location r = 2.5 in the simulations to 
be equivalent to 5 AU, then the time taken for a typical 
planetesimal to migrate a distance equal to its semimajor 
axis is t ~ 1 Myr. The stochastic torques were generated 
by a disk model that is approximately three times more 
massive than a minimum mass model. A lower mass disk 
would generate a correspondingly longer migration time. 
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It is clear that planetesimals in turbulent disks will 
undergo significant migration during the lifetime of the 
nebula. From the point of view of forming planets quickly, 
turbulence induced changes to the orbital elements can 
have both a positive and a negative effect. Increased 
mobility generally acts to favour forming planets more 
quickly. One of the major issues facing the theory of 
planet formation is forming the cores of the giant planets 
fast enough to accrete gas before dispersal of the nebula. 
Although estimates of the core formation time based on 
models of runaway growth suggest relatively short core 
formation times (e.g. Pollack et al 1996), simulations in- 
dicate that runaway growth does not continue all the way 
to the completion of core formation. Runaway growth 
slows down and enters a stage of 'oligarchic growth' when 
the eccentricities and inclinations of the planetesimals 
are pumped up by the larger embryos (Ida & Makino 
1993; Kokubo & Ida 1998). Simulations by Thommes et 
al (2003) indicate that core formation becomes difficult, 
in part due to gap formation in the planetesimal disk by 
the larger embryos and in part due to the dynamical exci- 
tation of the planetesimal disk, leading to cores of ~ few 
M ffl forming in a few million years. The isolation of cores 
may be alleviated by the increased mobility of embryos 
and planetesimals in a stochastic migration scenario. In 
addition, the existence of regions of the disk where the 
turbulence is weaker than in others may lead to the accu- 
mulation of planetesimals there that could speed up accu- 
mulation processes. 



The eccentricity and (and presumably) inclination 
driving by turbulence indicated by the right panel of fig- 
ure 4, however, may slow down the early stages of planet 
formation. In particular the runaway growth of planetes- 
imals will be affected, as the velocity dispersion of small 
bodies needs to be smaller than the escape velocity from 
the larger accreting objects. Indeed a large velocity dis- 
persion generated by turbulence may lead to destructive 
rather than accumulative collisions between planetesimals, 
leading to questions about how planetesimals form at all 
in a turbulent disk. It should also be noted that the large 
scale migration of icy planetesimals from the outer solar 
system into the inner solar system may lead to the inner 
planets being more enriched in volatiles than is observed. 



These problems may be alleviated if the stochastic 
torque values obtained in the MHD simulations performed 
here, using cylindrical disk models, turn out to be overes- 
timates when models including vertical stratification are 
computed. A vertical gradient may then occur in the rela- 
tive density fluctuations due to magnetic buoyancy effects 
(e.g. Stone et al 1996; Stone & Miller 2000), which reduce 
the torque fluctuations at the disk midplane. The effects of 
a vertical gradient in the ionisation fraction (e.g. Gammie 
1996) would act in a similar manner. 



9.2. Long term evolution of protoplanet orbits 

The situation regarding the long term evolution of plan- 
etesimals in turbulent disks is clear they will slowly dif- 
fuse throughout the disk. The situation regarding high 
mass, gap forming planets also seems to be clear. They 
will migrate inward at the effective viscous evolution of 
time of the disk, as in the standard type II migration pic- 
ture (Nelson & Papaloizou 2003, 2004). The situation re- 
garding low mass protoplanets remains ambiguous. The 
central question remains: can stochastic forces overcome 
type I migration and prevent at least some low mass pro- 
toplanets from falling into the star prior to accreting a 
gaseous envelope and becoming gas-giants ? 

The answer to this question depends on how effec- 
tive type I torques are in a turbulent disk, and whether 
the disk turbulence can generate fluctuations with the re- 
quired amplitudes and temporal behaviour to counterbal- 
ance the inward drift on planet formation time scales of 
~ 1 Myr. The simulations presented in sections 6.3 - 6.9 
and section 7.2 indicate that type I torques may be dimin- 
ished in turbulent disks compared with those that arise 
in laminar disks, but the evidence is far from conclusive. 
The power spectra presented in figure 13 show that the 
disk turbulence is capable of generating long term fluc- 
tuations with the appropriate amplitude to significantly 
affect type I migration, but there is no way of predicting 
whether such fluctuations can persist for time scales ap- 
propriate to planet formation itself. This question can only 
be answered by simulations whose run times are currently 
too long to be performed. 

The origin of the longer time scale fluctuations ob- 
served in the simulations remains unknown. We have sug- 
gested that global communication across the disk on the 
viscous time scale may feed into the fluctuations, con- 
tributing to their long term temporal evolution. The vis- 
cous communication time from the outer edge to the 
planet forming region in our model disks is ~ 8000 or- 
bits. For naturally occuring disks, that form through the 
collapse of molecular clouds and have radii in excess of 
100 AU, the viscous evolution time exceeds 1 Myr. It is 
not unreasonable to suppose that the nature of the turbu- 
lence in the planet forming region between 1-10 AU may 
be affected by global communication through the disk on 
this time scale. This may in turn feed into the torques 
experienced by forming protoplanets. 

Other effects that may arise in naturally occuring disks 
are the formation of surface density depressions or 'gaps'. 
Simulations show that variations in the magnetic stresses 
occur radially throughout the disk, and persistence of 
these can lead to the formation of 'gaps' (e.g. Hawley 
2000; Steinacker & Papaloizou 2002). Such features may 
play an important role in modifying the migration of low 
mass planets. The tendency of MHD turbulence to assist 
in gap formation for more massive planets has already 
been noted by Nelson & Papaloizou (2003) and Winters, 
Balbus & Hawley (2003). Radial variations in the strength 
of the disk turbulence may also play a role. The existence 
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of 'dead zones' has been suggested, particularly in the 
planet forming region around 5 AU (e.g. Gammie 1996; 
Fromang, Terquem & Balbus 2002), with an inner active 
zone maintained through thermal ionisation of alkalii met- 
als. A dominant outward flux of waves excited by a more 
active interior disk may provide a bias that helps stochas- 
tic torques counter balance inward type I migration. 

10. Conclusions 

In this paper we have presented simulations of low mass 
protoplanets embedded in turbulent, magnetised disks. 
The planet orbits are evolved, with forces due to gravita- 
tional interaction with the disk being included. The disk 
models are cylindrical, as the vertical component of grav- 
ity is neglected, but are otherwise designed to be similar to 
disks thought to be present around young stars during the 
planet formation epoch (i.e. disk thickness H/r = 0.07; 
disk viscosity a~5x 10~ 3 ). 

Planet masses, m p , ranging between zero and 30 M® 
were considered. The m p — 'planetesimals' were found 
to migrate through the disk in a manner similar to a ran- 
dom walk, and in some cases their eccentricity was found 
to increase to values e ~ 0.1. The random walk element 
of their evolution provides a means of maintaining signif- 
icant mobility of the planctcsimal swarm during planet 
formation. This may help reduce the problems of orbital 
repulsion and core isolation found during simulations of 
giant planet core formation (e.g. Thommes et al. 2003). 
However, the increase in eccentricity may present a prob- 
lem during the runaway growth stage of planet formation, 
which requires a dynamically cold disk of planetesimals 
(e.g. Wetherill & Stewart 1993). Further simulations that 
include gas drag and vertically stratified disk models are 
required to obtain a more accurate picture of this eccen- 
tricity growth. 

The protoplanets with finite masses that were studied 
also showed evidence of stochastic migration and eccen- 
tricity growth. The eccentricity evolution was such that 
the protoplanets with mass m p — 1 and 3 M® showed 
peak eccentricities of e ~ 0.08 and 0.045, respectively, 
whereas the m p = 10 and 30 M® planets had peak ec- 
centricities of e ~ 0.02. This is in broad agreement with 
expectations, as the interaction of low mass protoplan- 
ets with disks causes eccentricity damping due to interac- 
tion with material at coorbital Lindblad resonances (e.g. 
Artymowicz 1993). This damping is expected to scale lin- 
early with the protoplanet mass, leading to the observed 
behaviour. 

For all protoplanet masses in the range 1 < m p < 30 
M®, there was clear evidence of stochastic torques be- 
ing dominant over the type I migration expected to occur 
in laminar disks. We considered an ensemble of models, 
and while some planets showed a definite trend toward 
inward migration, there were examples for each mass con- 
sidered for which stochastic migration was the dominant 
effect. This behaviour arises in part because the turbulent 
fluctuations contain low frequency components whose am- 



plitudes are such that they can significantly modify type 
I migration over the simulation run times that are cur- 
rently feasible. In addition, we analysed the effect of su- 
perposing type I migration expected for laminar disks onto 
the stochastic migration obtained in turbulent disks. This 
analysis suggested that inward migration ought to be more 
apparent in the MHD simulations than was found to be 
the case, providing a hint that type I migration may be 
modified in turbulent disks. We note, however, that the 
simulations arc not conclusive on this latter point. 

A number of important issues remain to be addressed. 
These include a determination of the role of gas drag in 
modulating the effects of eccentricity growth of planetes- 
imals, and the effects of vertical disk structure in influ- 
encing the strength of the stochastic torques at the disk 
midplanc. The question of whether stochastic torques can 
prevent type I migration of protoplanets over planet for- 
mation time scales still remains to be settled, and can 
be partially answered by performing longer simulations. 
The potential modification of type I migration in turbu- 
lent disks can be addressed by customised simulations that 
examine the exchange of angular momentum between disk 
and planet These will be the subjects of future publica- 
tions. 
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